Interpreting the Atmospheric Composition of Exoplanets: Sensitivity to Planet Formation Assumptions

Constraining a planet's formation location based on its elemental abundance ratios (e.g., C/O) requires that the elemental composition of the protoplanetary disk is known. A good starting point may be to assume that the protoplanetary disk has an elemental composition that is identical to that of the host star (planet formation may deplete stellar photospheres in metals with respect to the disk, however; see Chambers 2010; Bitsch et al. 2018a; Booth & Owen 2020). It is therefore crucial to have knowledge about the host star's abundances that is as complete as possible, in excess of just [Fe/H], or to at least use existing scaling relations to approximate the stellar metal content for elements other than iron (e.g., Bitsch & Battistini 2020). The disk composition and assumed mass also set an upper limit on the amount of metals that a planet may accrete during its formation (e.g., Baraffe et al. 2008). Alternatively, the retrieved metal content of an exoplanet may also be used to place a lower limit on the disk metal content, or even total mass, analogous to the concept of the minimum-mass solar nebula (e.g., Hayashi 1981).

Moreover, the disk's physical and thermal structure is important to set the radial fractionation of elements into different molecular species, in solid and gaseous form, that can be accreted by a forming planet. The disk structure depends on the assumed (dust) opacities and therefore the dust evolution (Schmitt et al. 1997; Birnstiel et al. 2016; Savvidou et al. 2020). Moreover, disks are viscously (Lynden-Bell & Pringle 1974) evolving (or due to photoevaporation and disk winds; see Clarke et al. 2001; Bai et al. 2016; Suzuki et al. 2016; Chambers 2019), prone to different instabilities (e.g., Flock et al. 2017; Klahr et al. 2018), and will be affected by the presence of (especially massive) planets that may induce spiral density perturbations, lead to the formation of vortices, or open gaps (e.g., Lin & Papaloizou 1986; Crida et al. 2006; Lobo Gomes et al. 2015; Pinilla et al. 2015; Binkert et al. 2021).

Another important effect determining the disk gas composition is the evaporation of pebbles inside of ice lines, which may significantly increase the local volatile content of gas in the disk (Piso et al. 2015; Booth & Ilee 2019; Schneider & Bitsch 2021a). This could lead to planets being enriched in volatiles much more than expected from pure gas accretion in the classical Öberg et al. (2011) setup. It is important to mention that the dynamics of pebbles is likewise determined by the disk structure. The structure sets the pebbles' growth rates, their Stokes parameters (i.e., drift speeds), and may trap pebbles into local pressure maxima (e.g., Paardekooper & Mellema 2006; Ormel & Klahr 2010; Birnstiel et al. 2012; Lambrechts et al. 2014).

The chemical composition of the protoplanetary disk (e.g., Henning & Semenov 2013) is of central importance for determining the composition of planetary building blocks, that is, the disk's gas and solid phases. In Öberg et al. (2011) and the more recent Öberg & Wordsworth (2019) study, a simplified and static disk chemical model is assumed. In practice, the disk's chemical composition will evolve because both the disk gas and the volatile ices on grain surfaces will undergo chemical processing (e.g., Eistrup et al. 2016, 2018; Molyarova et al. 2017). This means that chemical reactions between atoms and molecules in both gas and ice may alter which molecular species are the dominant carriers of elements such as C, O, and N over time. This is an important effect and is expected to alter the inferred formation history of a planet, as has been pointed out by Eistrup et al. (2016). Examples are the conversion of CO into CO₂ ice over time, or the conversion of N₂ gas into NH₃ ice (Semenov & Wiebe 2011; Molyarova et al. 2017; Eistrup et al. 2018). Another example of the importance of disk chemistry on planet formation are the processes that lead to the observed carbon depletion in the inner solar system (e.g., Mordasini et al. 2016; Cridland et al. 2019), which may be attributable to the irreversible chemical destruction of carbon grains within a disk's so-called soot line, or connected to chondrule formation (e.g., Kress et al. 2010; Gail & Trieloff 2017; van' t Hoff et al. 2020; Li et al. 2021).

The disk chemistry itself is sensitive to many processes, such as stellar evolution (Miley et al. 2021), or the cosmic-ray ionization rate (Eistrup et al. 2016; Schwarz et al. 2019). Moreover, whether or not the initial composition of the disk matter is molecular, that is, "inherited" from the composition of the natal molecular cloud, or elemental ("reset" scenario) can strongly influence the C/O values (Eistrup et al. 2016). The disk's physical structure may also have an impact on its chemical evolution. As an example, we highlight the effect of the self-shadowing of the disk, allowing for nominally too cool compositions to occur closer to the star than otherwise expected (Ohno & Ueda 2021).

The idea that planet formation may be constrained through planet composition, as presented in Öberg et al. (2011), conceptually boils down to comparing the planetary C/O value and total metal enrichment to the C/O of the disks' solid and gaseous phases as a function of orbital distance from the star. However, planet formation involves and connects many complex processes. This means that a forming planet cannot accrete arbitrary amounts of gas or solids at (or from) arbitrary positions in the disk.

For example, a limiting factor for planets forming via the core accretion paradigm (e.g., Mizuno 1980; Pollack et al. 1996), specifically when accreting pebbles, is the pebble isolation mass. Pebble accretion, which is the accretion of roughly centimeter-sized solids by the forming planet (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012), can only dominate the solid accretion process until the growing planet reaches this isolation mass M_iso (e.g., Lambrechts et al. 2014; Ataiee et al. 2018; Bitsch et al. 2018b), after which pebble accretion stops. This is because the planet induces the formation of a pressure bump in the disk, exterior to its orbit, which traps the inward-drifting pebbles. The isolation mass places an upper limit on the refractory content of a planet. A refractory content higher than allowed by this concept of M_iso could point to the importance of accreting planetesimals (e.g., Mordasini et al. 2016; Brügger et al. 2020), unless the planet formed very close to its star, within the refractories' ice lines (Schneider & Bitsch 2021b).

For planets growing in situ, a planetesimal isolation mass was established by Lissauer & Stewart (1993). It is caused by the planet depleting the local reservoir of planetesimals within the zone of its gravitational influence (the so-called "feeding zone"). In contrast to the pebble isolation mass, the planetesimal feeding zone increases with planet mass. The different ways of how a planet's refractory content and total mass scale in planetesimal and pebble accretion models might thus allow us to put limits on the contributions of the different solid reservoirs for a given planet. The planetesimal isolation mass can also be overcome via (giant) impacts or if a protoplanet migrates into regions of the disk still containing planetesimals (Alibert et al. 2005).

Moreover, the accretion of gas by the growing planet is a three-dimensional process (e.g., D'Angelo et al. 2003; Ayliffe & Bate 2009; Szulágyi et al. 2014; Ormel et al. 2015; Schulik et al. 2019, 2020). This may be especially important, as the gas composition and C/O value are thought to be changing above the midplane of the disk (e.g., Molyarova et al. 2017; Cridland et al. 2020a). Similar to the solid accretion processes, the amount of gas a planet can accrete during formation is limited. In contrast to the solids, however, the ultimately limiting factor is the lifetime of the protoplanetary disk. Once a planet enters runaway accretion, it will accrete gas as quickly as can be provided by the viscously evolving disk, until the disk dissipates. More specifically, it has been shown that gas delivery to the planet can be severely limited by gap formation, which in turn is controlled by the disk's viscous resupply of the planetary feeding zone (e.g., Lubow et al. 1999; Ayliffe & Bate 2009; Lissauer et al. 2009; Bergez-Casalou et al. 2020; Schulik et al. 2020). For embedded planets below the gap opening mass, gas accretion may be limited if radiative cooling is counteracted by hot inflowing gas from the ambient disk. In this case even gas that enters the planetary Hill sphere is not accreted (so-called "recycling"; Cimerman et al. 2017).

Another important process of planet formation is migration. The migration history of a planet is critical to its final composition because it determines where in the disk it accretes material. Because the Öberg et al. (2011) approach strives to ultimately constrain formation locations with respect to the disk ice lines, migration may be less of a problem if a forming planet did not migrate across ice lines. Interestingly, planets forming in the inner disk may actually be trapped at locations just beyond the water ice line (e.g., Bitsch & Johansen 2016; Cridland et al. 2016; Müller et al. 2021). There also exist other traps not connected to ice lines, such as the disk's dead-zone inner edge (Bitsch et al. 2014; Cridland et al. 2016). Qualitatively speaking, and when not currently trapped, planets are expected to migrate either by fast type I migration or via slower type II migration, the latter of which ensues once the planet is massive enough to open a gap in the disk (see, e.g., the review by Baruteau et al. 2016). Quantitatively, there is an ongoing debate about the actual magnitude of the torques (and therefore speed of migration), for example, for type II migration (e.g., Dürmann & Kley 2015; Robert et al. 2018).

Further complicating the picture are N-body interactions between the forming planets. This process is now regularly included in models of planet formation and population syntheses but comes at an increased numerical cost (e.g., Alibert et al. 2013; Chambers 2016; Lambrechts et al. 2019; Emsenhuber et al. 2021; Izidoro et al. 2021). While increasing the complexity of describing planet formation, N-body interactions may represent an alternative avenue for producing hot Jupiters, which may have been scattered in by planets farther out (e.g., Jurić & Tremaine 2008; Bitsch et al. 2020). N-body interactions can also alter the composition of a planet via giant impacts. The amount of solids brought into a giant planet via such impacts might be substantial or even dominant compared to the amount accreted from small bodies like planetesimals or pebbles (e.g., Emsenhuber et al. 2021; Ginzburg & Chiang 2020; Ogihara et al. 2021). If the impactors originally formed in clearly different regions of the disk than the forming planet, this would blur the meaning of a well-defined formation location of a planet. Accounting for the effects of N-body interactions when trying to invert planet formation thus seems challenging. Interestingly, it has recently been shown that while N-body interactions tend to randomize the process of planet formation, machine-learning techniques such as random forests still allow us to predict the outcome of planet formation quite accurately (Schlecker et al. 2021). In this work the authors show that the initial parameters of the planet formation model described in Emsenhuber et al. (2021), mainly the initial location of the planetary embryo and the dust mass of the disk, may be used to predict which class a forming planet will belong to (super-Earths, Neptune-like, giant planets, etc.). These classes also correspond to certain orbital distances and planetary compositions.

The discussion here focused mostly on planet formation via the core accretion paradigm. Other aspects are of importance if a planet forms via gravitational instability (GI; see Boss 1997). The disk structure (and thus formation environment of the planet) will be quite different in this case. This is because GI planets may form early, when the disk is still massive, in the outer parts of the disk (Boss 2021; Schib et al. 2021). We note that while GI is classically regarded as a way to produce wide-separation gas giant planets, it has also been suggested to allow for the formation of less massive, small-separation planets (Nayakshin 2010), caused by fast inward migration after formation and the associated mass loss, dubbed "tidal downsizing."

When aiming to constrain planet formation based on the results of atmospheric abundance characterizations, one must make assumptions on how the atmospheric composition relates to the bulk composition of a planet. It has been pointed out that planets growing via core accretion may have a layer of heavily metal-enriched gas above their solid cores, due to the evaporation of pebbles and planetesimals that are destroyed when entering the hot planet's proto-atmosphere (e.g., Mordasini et al. 2016; Helled & Stevenson 2017; Brouwers et al. 2018; Brouwers & Ormel 2020). The recent constraints by the Juno spacecraft on the interior of Jupiter are consistent with this assessment, pointing to the existence of a dilute core (Wahl et al. 2017; Debras & Chabrier 2019). We note that the exact distribution of metals in Jupiter' interior is difficult to explain, however, and may involve a formation process stretching over 2 Myr (see Helled et al. 2022, for a recent review). What is more, the planet bulk metallicity constraints for transiting planets derived in Miller & Fortney (2011), Thorngren et al. (2016), and Thorngren & Fortney (2019) tend to be higher than the metallicities reported for planetary atmospheres (albeit with large uncertainties; see Welbanks et al. 2019), making an enrichment of the interior with respect to the atmosphere a likely scenario, and allowing for a first estimate regarding the efficiency of planetary mixing.

The question that naturally arises from these findings is how representative the inferred atmospheric composition is of the bulk of the planet, even for giant planets, where the gas dominates the mass budget. For this it needs to be understood if and how well the metals can be mixed throughout the planetary envelope. Whether this happens at all is not clear, as a gradient in metallicity (and therefore mean molecular weight) may stabilize the planetary interior against convective motions (e.g., Ledoux 1947; Chabrier & Baraffe 2007; Leconte & Chabrier 2012). How a solid core or compositional gradients tend to mix throughout the planetary envelope for gas giant planets, often taking Jupiter as an example, is currently being investigated (Vazan et al. 2015, 2016; Moll et al. 2017; Vazan et al. 2018; Müller et al. 2020; Ormel et al. 2021). The results reported in these studies do not yet agree, predicting either fully mixed envelopes or ones where a metal gradient (and core) persists in the planet. In addition to imperfect mixing, rainout processes also likely play an important role, potentially depleting both solar system and exoplanet atmospheres in metals (e.g., Spiegel et al. 2009; Wilson & Militzer 2010).

From the discussion above it becomes clear that the metal enrichment inferred from the atmospheric retrievals, which serve as an input for any formation analysis, is only a lower limit for the true planetary metal enrichment. As long as the metals are locked into the invisible interior of the planet homogeneously, and not selectively, this may still allow us to constrain the solids' location of origin in the protoplanetary disk, as long as the atmosphere is still enriched enough for solid accretion to be the dominating factor. In the case where the relative elemental composition in the atmosphere (except for H and He) is different from the deep interior, or when the atmosphere is depleted with respect to the deeper interior to the point where it mimics the atmospheric metallicities expected from pure gas accretion, this poses a problem. Depending on how planet formation ensued, the former case may still trace the origin of the solids that were accreted toward the end of the formation process. The latter case could be resolved by checking the relative atomic abundances in the presumed metal-poor planet. If refractory atomic species are relatively abundant, this may point to a dominant accretion of solids that are mostly locked into the planet's interior. An interesting recent discussion of how to constrain planet formation in the case where the atmospheric and planetary bulk compositions differ can be found in Helled et al. (2022), who come to similar conclusions. We note again that these avenues for analyzing the origin of the metals in a planet's atmosphere will be further complicated if volatile-rich gas from evaporated pebbles was accreted by the planet.

The atmospheric composition of an exoplanet may also evolve owing to atmospheric evaporation, or by secular enrichment of infalling comets and asteroids. Evaporation is especially important for close-in low-mass planets (e.g., Jin & Mordasini 2018). The atmosphere may be partially or fully lost owing to thermal or nonthermal processes, where the thermal escape separates into the regimes of Jeans escape or hydrodynamical escape, depending on the local thermal state of the atmosphere (see, e.g., Barman 2018). For the atmosphere to become relatively enriched in metals by evaporation, two criteria have to be met. First, the atmosphere needs to be of low enough mass to allow for a significant amount to be lost. Second, the atmospheric escape process needs to preferentially retain the heavier atmospheric species, for which the atmosphere needs to be in the Jeans escape regime (e.g., Bourrier & Lecavelier des Etangs 2018). In the hydrodynamic escape regime the heavier metal species would be lost together with hydrogen. We note that this transition is gradual, however, and that there can be mass fractionation in hydrodynamic outflows as well, depending on the magnitude of the total mass flux (e.g., Hu et al. 2015). These authors also report that such outflows can selectively deplete the atmospheres of Neptune- and sub-Neptune-mass planets in hydrogen over multiple Gyr, provided that the initial atmospheric mass is small enough (<10⁻³ of planetary mass). For gas giant planets evaporation may thus be a less relevant process for changing the atmospheric composition. An extreme case that is worth mentioning is the Roche lobe overflow that may affect the closest-in gas giant planets. This process could strip away the upper gaseous envelope, potentially revealing the more metal-enriched layers below. Roche lobe overflow has been discussed as the potential origin of LTT 9779b, a planet in the hot Neptune desert (Jenkins et al. 2020). By extension, significant atmospheric evaporation may lead to similar outcomes if a compositional gradient is present in the atmosphere.

Another possibility of atmospheric evolution is the secular contamination of the planetary atmospheres by infalling comets or asteroids. Here the frequency of cometary impacts and the persistence of the enrichment they cause in the atmosphere need to be estimated. Turrini et al. (2015) find that the additional water a comet may deposit in the visible atmosphere of the hot Jupiter HD 189733b would have to persist for 500–5000 yr before being removed, assuming impacts of kilometer-sized comets every 20–200 yr; otherwise, no significant enrichment is possible. A quantitative assessment of cometary enrichment requires an estimate of the persistence timescale, however. As this appear to be lacking from the literature, we present a simple first-order analysis below. We start by approximating local mixing timescales as ${\tau }_{\mathrm{mix}}={H}_{P}^{2}/{K}_{{zz}}$, where H_P is the planetary pressure scale height and K_zz the atmospheric eddy diffusion coefficient, and find

$$\begin{eqnarray}&&{\tau }_{\mathrm{mix}}=0.3{\left(\displaystyle \frac{{H}_{P}}{300\ \mathrm{km}}\right)}^{2}\left(\displaystyle \frac{{10}^{8}\ {\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}}{{K}_{{zz}}}\right)\ \mathrm{yr}.\end{eqnarray} \tag{ 1 }$$

A K_zz value of 10⁸ cm² s⁻¹ within the radiative zone is well within the estimates by passive tracers reported from model calculations for HD 189733b and HD 209458b (Parmentier et al. 2013; Agúndez et al. 2014). For self-luminous planets the correct value to be chosen is unclear, but K_zz = 10⁵ cm² s⁻¹ may at least be a useful lower limit (Ackerman & Marley 2001). In the deeper regions of the atmosphere convective overshoot may lead to higher values for K_zz , smoothly transitioning toward the value expected for fully convective atmospheres as the radiative–convective boundary (RCB) is approached (e.g., Ludwig et al. 2002; Helling et al. 2008).

The τ_mix estimate of less than a year in Equation (1) thus shows that mixing potentially proceeds on faster timescales than those quoted by Turrini et al. (2015), at least locally, meaning that any material added by cometary impacts should be mixed away quickly. This does not preclude a slower, homogeneous enrichment of the radiative atmosphere by cometary impacts over time. However, the question is how quickly this enrichment will be removed into the bulk interior of the planet by entrainment into overshooting convective blobs at the RCB. By extension, the enrichment of the visible atmosphere may be lower if the impactors deposit their metals below the RCB: for the bulk of Jupiter, for example, the mixing timescale has been estimated to be at most a few years (Debras & Chabrier 2019).

Assuming diffusive mixing in a 1D atmosphere, we derive the following expression for the increase in mass fraction X of a certain species in the planetary atmosphere, due to cometary impacts:

$$\begin{eqnarray}&&{\rm{\Delta }}X=\frac{\dot{M}g}{4\pi {R}_{{\rm{P}}}^{2}}\frac{{H}_{P}^{2}}{{K}_{{zz}}}\frac{{P}_{\mathrm{RCB}}-{P}_{{\rm{i}}}}{{P}_{\mathrm{RCB}}{P}_{{\rm{i}}}}.\end{eqnarray} \tag{ 2 }$$

Here $\dot{M}$ is the mass accretion rate of comets, g the planetary gravity, R_P the planetary radius, H _P the atmospheric pressure scale height, and P_i the average location of the destruction of impacting comets in units of pressure. P_RCB is the location of the RCB or, alternatively, the location where the deep mixing becomes fast enough to make the local mass fraction equilibrate to the average of the planetary interior (which is assumed to be well mixed). In any case it must hold that P_i < P_RCB for this expression. A derivation can be found in Appendix A. When assuming pure-water comets and tracking the change in the mass fraction of water, a relative enrichment of ΔX/X₀ = 3 × 10⁻⁴ to 3 × 10⁻² is found, where X₀ is the planet's bulk water mass fraction. See Appendix A for more details and which values to assume for the various quantities. We therefore conclude that secular cometary enrichment of the planetary atmosphere may be unlikely for giant planets, but a better modeling of the process, for example, the lower boundary treatment when using a 1D diffusion approximation, is needed.

Formation models produce synthetic planets, or populations thereof, starting from a physical model and set of formation parameters. These can be, for example, the initial disk mass, the disk composition, and the starting position of a planetary embryo in a disk. When attempting to constrain planet formation based on measured planetary compositions, the formation models need to be inverted, because planet composition is an outcome of the formation models. More specifically, if assuming that planetary compositions

$$\begin{eqnarray}&&{ \mathcal C }=({\rm{O}}/{\rm{H}},{\rm{C}}/{\rm{H}},{\rm{N}}/{\rm{H}},\mathrm{Fe}/{\rm{H}},\mathrm{Si}/{\rm{H}},{\rm{P}}/{\rm{H}},{\rm{S}}/{\rm{H}},...)\end{eqnarray} \tag{ 3 }$$

can be measured with a given uncertainty, we are then interested in the probability distribution of formation parameters ϑ of a formation model ${ \mathcal M }$, given this measurement. Using Bayes's theorem, this can be written as

$$\begin{eqnarray}&&P(\vartheta | { \mathcal C },{ \mathcal M })=\displaystyle \frac{P(\vartheta | { \mathcal M })P({ \mathcal C }| \vartheta ,{ \mathcal M })}{P({ \mathcal C }| { \mathcal M })}.\end{eqnarray} \tag{ 4 }$$

Here $P(\vartheta | { \mathcal M })$ is the prior probability of ϑ before considering any data ${ \mathcal C }$, while $P({ \mathcal C }| \vartheta ,{ \mathcal M })$ is the likelihood for observing ${ \mathcal C }$, given that ϑ is true. In practice, this may be written as, for example,

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}P({ \mathcal C }| \vartheta ,{ \mathcal M }) & = & -\displaystyle \sum _{i=1}^{{N}_{\mathrm{species}}}{\left[{{ \mathcal C }}_{i}-{{ \mathcal C }}_{i}^{{ \mathcal M }}(\vartheta )\right]}^{2}/2{\rm{\Delta }}{{ \mathcal C }}_{i}^{2}+\mathrm{cst}\\ & = & -{\chi }^{2}/2+\mathrm{cst},\end{array}\end{eqnarray} \tag{ 5 }$$

where N_species is the number of measured atmospheric species, ${{ \mathcal C }}_{i}^{{ \mathcal M }}(\vartheta )$ is the formation model prediction of the planet abundance of species i, and ${\rm{\Delta }}{{ \mathcal C }}_{i}$ are the measurement uncertainties. Here we choose a simple form of the likelihood for clarity, assuming that the measured abundances of different atmospheric constituents are independent and follow a Gauss distribution. In general, the functional form of the log-likelihood can be arbitrarily complicated. For example, the abundance posterior of an atmospheric retrieval may be used directly, which can be approximated by, say, a Gaussian mixture model. In our application in Section 3 we chose an intermediate step, accounting for the covariance between the atmospheric oxygen and carbon content. We note that such an inversion process does not necessarily only need to consider elemental abundances as input measurements. Any observed property of a planet, such as its orbital parameters, could in principle be included in this analysis, as long as it is predicted by a formation model.

In practice, we will compute samples of the target distribution $P(\vartheta | { \mathcal C },{ \mathcal M })$ by numerically integrating the numerator of the right-hand side of Equation (4) using the so-called nested sampling method (Skilling 2004). In short, nested sampling is a Monte Carlo technique to integrate functions in highly dimensional model parameter spaces. Here it is the space spanned by the formation parameters. When integrating the numerator of Equation (4), the integral value resulting in the so-called model evidence, nested sampling will automatically generate samples of our target distribution. In principle, model evidence ratios allow us to distinguish between different formation models. However, as long as we cannot invert full state-of-the-art formation models, this may be possible only when considering sets of assumptions that lead to wildly different outcomes, and where one clearly represents a better fit. We use the PyMultiNest (Buchner et al. 2014) package for inverting formation models, which is a Python wrapper of the MultiNest code (Feroz et al. 2019). A schematic illustrating our approach is shown in Figure 2.

Zoom In Zoom Out Reset image size

3.1.1. Inverting the Öberg et al. (2011) Model

As a first toy model we will use arguably one of the simplest formation models conceivable when aiming at studying the usefulness of atmospheric C/O constraints. For this we follow the setup presented in Öberg et al. (2011) and Öberg & Wordsworth (2019), assuming their static protoplanetary disk model as described for the young solar nebula, but with the ice lines and composition adapted to the planet system of interest. For every volatile species we define a constant mass fraction in relation to the total disk mass. Inside of its ice line the volatile species is in the gas phase, whereas outside it is in the solid phase. This leads to the well-known step-like behavior of the C/O values in the solid and gaseous phase of the disk, as shown in Figure 1. More details on our implementation can be found in Appendix B.1.

We then assume that the planet formation process can be fully described by a set of four parameters, which ultimately allow us to map to the bulk composition of the planet:

$$\begin{eqnarray}&&\vartheta =({{\boldsymbol{M}}}_{{\rm{P}}},{M}_{\mathrm{solid}},{a}_{\mathrm{solid}},{a}_{\mathrm{gas}}),\end{eqnarray} \tag{ 6 }$$

where M_P is the total planet mass and M_solid is the mass of the solids (refractory species and volatile ices) accreted by the planet. It holds that M_P = M_gas + M_solid, which we use to determine the amount of gas a planet accreted. The parameters a_solid and a_gas denote the orbital distances where the solids and gas were accreted, respectively.

For a given value of ϑ our setup then uses the disk model to determine the planet's accretion location with respect to the disk ice lines. Species in the gas phase will be used to determine the composition of the accreted gas; the analogous is done for the solids. Together with the mass fractionation between solids and gas in the planet, this determines the planet's composition ${ \mathcal C }$.

3.1.2. Adding Chemical Time Evolution

The Öberg et al. (2011) disk setup is convenient for conceptually studying the usefulness of C/O. As discussed in Section 2, there are many complicating factors that make a true inversion from C/O to planet formation parameters extremely challenging. In our second toy model setup we single out one of these processes and study how chemical evolution in the protoplanetary disk changes inferences when compared to the static disk model.

As discussed before, chemical reactions may process the initial disk abundances over time, shifting carbon and oxygen atoms to different chemical species. For example, CO, which is very volatile, can be removed from the gas by surface reactions, processing it into less volatile CO₂ (Molyarova et al. 2017; Bosman et al. 2018; Eistrup et al. 2018; Schwarz et al. 2018). If a large fraction of CO gas in a given region of the disk midplane undergoes such processing into CO₂ ice, then an amount of elemental C and O equivalent to that initially carried in CO will be processed from the gas into the ice in this region. This is one example of how chemical processing in disks can alter the elemental partitioning of, for example, C and O, in the gas and ice. These processes are thus expected to influence inferences on the location of planet formation as first discussed in Eistrup et al. (2018).

To study the effect of disk chemical evolution in formation inversions, we replaced the static disk abundance model with a time-dependent model. For this we calculated the evolution of the disk chemical composition using the ANDES code, which describes a quasi-stationary 2D axisymmetric protoplanetary disk (Akimkin et al. 2013; Molyarova et al. 2017). ANDES solves for the time-dependent chemical composition of the disk with a detailed description of grain surface and gas-phase reactions; see Appendix B.2 for more details. For the initial disk setup we use the abundances given in Table 2, that is, equal to the disk abundances used for our static disk model. Although the chemical model includes other elements, we only consider H-, He-, C-, O-, and N-bearing species in the calculations. The refractories are considered to be chemically inert and condensed at all times. An example for the resulting time evolution of the C/O ratio in the disk is shown as the yellow to dark-blue lines in Figure 3. The process of planet formation is then modeled in the same way as in the static disk case, but with the difference that the composition of the accreted gas and solids is taken from the disk's chemical evolution calculations. Formally this also allows us to add the formation time as an additional parameter to be constrained, but here we chose to initially only study and compare the inversion outcomes for different times during chemical evolution.

Zoom In Zoom Out Reset image size

3.1.3. Studying the Effect of Pebble Accretion

As discussed in Section 1, deriving accurate and precise C/O ratios from current exoplanet observations is challenging, but this will likely change with JWST. Using high-resolution spectrographs, or the interferometric GRAVITY instrument at the Very Large Telescope, the community is already starting to derive precise C/O values using ground-based instruments. This has to do with the data but also with the use of state-of-the-art retrieval techniques (Brogi & Line 2019; Gandhi et al. 2019; Mollière et al. 2020; Line et al. 2021; Pelletier et al. 2021). Below we will make use of the atmospheric composition derived for HR 8799e from GRAVITY observations (Mollière et al. 2020) and try to constrain how its derived formation history changes when disk chemical time evolution or pebbles are included.

In order to run the formation inversion for HR 8799e, a disk elemental composition needs to be assumed for HR 8799. HR 8799 is a λ Boötis–type star; this means that the abundances measured for its iron-peak elements are subsolar, with values of [Fe/H] = −0.55 ± 0.10 (Sadakane 2006) or [Fe/H] = −0.52 ± 0.08 (Wang et al. 2020) having been inferred for iron specifically. A similar depletion is expected for Mg, Si, and other massive iron-peak elements, while the abundances of elements typically found in volatile species (C, N, O) are expected to be close to solar (e.g., Paunzen 2004). Indeed, the latest analysis of Wang et al. (2020) inferred [C/H] = 0.11 ± 0.12, [O/H] = 0.12 ± 0.14, and (C/O)/(C/O)_⊙ = 0.96 ± 0.19, which are all consistent with solar but slightly enriched in C and O. According to Wang et al. (2020), the most likely explanation for the observed composition of HR 8799 is recent accretion of volatile-rich material onto the outer layers of the star, for example, from an evaporating hot Jupiter, or of volatile-rich ices scattered into the inner system by the four HR 8799 planets. We will therefore use the composition of the solar system from Öberg & Wordsworth (2019) as our nominal abundance model because it could be unlikely that the star has a bulk elemental composition identical to its observed photospheric values. When relevant, we will also report on how our results change if taking the λ Boo abundances of HR 8799 at face value, however. The ice line positions in the HR 8799 disk are set to the values derived from the ANDES chemical evolution model, at t = 0.

For HR 8799e we make use of the atmospheric retrieval results reported in Mollière et al. (2020). Of relevance for the formation inversion are the derived values for the planetary mass, as well as the atmospheric metallicity and C/O ratio. As the mass is spectroscopically determined, it has large error bars. However, it still results in a constraint on the total amount of solids that the planet incorporated, which can also be estimated from multiplying the atmospheric metallicity by the planetary mass. More specifically, multiplying the planet mass by the inferred atmospheric metallicity results in a lower limit on the metal mass (in solid or gaseous state) that a planet accreted, which may be dominated by solids in cases of high metallicities. We used the actual posterior distributions on planetary gravity and radius from the spectral retrievals to construct the inversion prior for the planetary mass (effectively corresponding to a 1σ upper limit of 14 M_Jup). We also study the effect of using a tighter mass constraint in the pebble inversion scenario.

Converting the atmospheric C/O ratio for use in the inversion method requires special care. In the spectral retrievals the metallicity is used as a free parameter to scale all elemental abundances except H and He, after which the C/O ratio is set by scaling O. In the formation model the formation location and relative gas-to-solid accretion fraction set the O, C, and refractory metal content, from which a C/O ratio can be calculated. Therefore, the C content is no longer strictly coupled to the refractory metal content, in contrast to the spectral retrievals that we use as input for our formation retrievals. This inconsistency has to be kept in mind when we use the C/O and metallicity of the spectral retrieval to obtain atmospheric C/H and O/H values and compare these to the C/H and O/H predicted by the formation model. In general, independently constraining C/H and O/H in atmospheric retrievals is the better avenue for running formation inversion studies. We note that it is also important to account for the amount of oxygen that has been sequestered into atmospheric clouds. Because the C/O constraints from Mollière et al. (2020) include this effect, the atmospheric retrieval results can be used without modification. For C/O constraints from retrievals that constrain absorber abundances independently, corrections would need to be applied.

Moreover, sampling the C/O and metallicity posterior of the spectral retrieval leads to tightly correlated C/H and O/H values, and we take this into account by fitting it with a two-dimensional Gaussian distribution. The distribution's covariance matrix is then used to describe the uncertainties of C/H and O/H during the inversion process. If this were not done, the independent uncertainties in C/H and O/H (as obtained from the diagonals of the covariance matrix) would allow for a spread in C/O values much larger than obtained from the spectral retrieval, rendering a formation inversion meaningless. Taking the spectral retrieval results from Mollière et al. (2020), we used [Fe/H] = 0.48 ± 0.27 and C/O = 0.60 ± 0.075. The resulting 2D distribution of C/H and O/H is shown in Figure 4.

Zoom In Zoom Out Reset image size

Lastly, because our adapted Öberg & Wordsworth (2019) setup for the disk abundance led to a slightly subsolar C/O ratio (0.52 instead of 0.55), we scaled the planetary value of 0.6 by 0.52/0.55 prior to generating the planetary C/H, O/H distribution as input for the formation inversion process. This is done to conserve the relative distance in C/O, with respect to solar abundances, of the input planetary composition.

3.2.1. Static Disk Chemistry

We start the formation inversion for HR 8799e using the disk model à la Öberg et al. (2011) and Öberg & Wordsworth (2019). The magenta ellipse in the top left panel of Figure 4 shows the distribution of sampled C/H, O/H pairs of the formation inversion, indicating a good fit. The most interesting result from the formation inversion is shown in the left panel of Figure 5: here we plot the 1D posterior of a_solid, which is the location where the solids that are enriching the planet originated or were accreted.

Zoom In Zoom Out Reset image size

The inversion process finds a clear preference for the solids to stem from outside the CO ice line, or from within the H₂O ice line. This is intuitively easy to understand: because the spectral retrieval resulted in a superstellar atmospheric metallicity and a C/O ratio consistent with stellar, this means within the Öberg et al. (2011) model that the planet must have accreted solids of stellar C/O ratio. The metal content of the gas, being stellar or substellar, is not high enough to offset the C/O ratio set by the accretion of solids. Accreting solids of stellar C/O is possible at orbital distances outside of the ice lines of all major carbon- and oxygen-carrying species, the outermost being CO (see Figure 1). For HR 8799e's current orbital position, ∼15 au (Wang et al. 2018), this could mean that the planet underwent some orbital migration after solid accretion, as the CO ice line for a young A5 host star such as HR 8799 is expected to be around 35 au. To illustrate this, we overplot today's orbital location of HR 8799e in the left panel of Figure 5.

Alternatively, due to a high organic carbon content of the refractories in our toy model, a roughly stellar C/O ratio is also attainable if the planet formed within the water ice line and then migrated or scattered outward to its current orbital position. This formation channel for distant giant planets was proposed, for example, in Marleau et al. (2019). While this is an intriguing result, we stress that it is dependent on the disk compositional model we assume and the formation model used in general. We also note that the high carbon content of the refractories in the inner disk may be unlikely; see our discussion in Section 2. The 1D and 2D projection of the full posterior of the formation inversion is shown in the left panel of Figure 12 and discussed in Appendix C.

We also carried out inversions using the λ Boo composition of the star for the disk. This was done by increasing the oxygen and carbon abundance by 30% and decreasing the iron and silicate content of the refractory material by 70%. The oxygen no longer bound in silicates was added to H₂O, which is the dominant reservoir of oxygen in the protoplanetary disk. CO, the second most abundant oxygen reservoir, should not change because the carbon content of the disk is not changed when applying the depletion of the iron-peak elements. Finally, the oxygen abundance is increased until C/O = 0.54 is reached, which is the value reported in Wang et al. (2020). We find that the most likely location of the origin of the accreted solids is still outside the CO ice line. The formerly second likely location, inside the H₂O ice line, vanished: the solid C/O there is exclusively set by the refractory species, which have a much higher C/O value of 2.6 now because of the strong silicate depletion. The associated posterior is shown in the right panel of Figure 12, in Appendix C.

3.2.2. Chemical Disk Evolution

3.2.3. Pebble Drift and Evaporation

In this section we model the effect of pebble drift and evaporation on the formation inversion of HR 8799e. This process can be crucial for setting the composition of forming planets. When neglecting pebble drift, planets whose atmospheric metal content is set by gas accretion are generally expected to have substellar metallicities and superstellar C/O values. In contrast, planets with an atmospheric metal enrichment dominated by solid accretion may have superstellar metallicities but substellar C/O ratios (e.g., Öberg et al. 2011; Madhusudhan et al. 2014, 2017; Mordasini et al. 2016). In the case of pebble drift, however, evaporation of pebbles inside of the CO, CO₂, and potentially the CH₄ ice lines can lead to disk gas that is significantly enriched in these species, allowing for superstellar metallicities and C/O ratios in the disk's gas phase and therefore in the atmospheres of planets (e.g., Booth et al. 2017; Schneider & Bitsch 2021a, 2021b).

For setting up the chemcomp pebble disk model, we used the same initial disk surface density and temperature structure as for the disk's chemical evolution case. Likewise, the initial disk composition was fixed to the one described in Table 2. For the disk viscosity we chose an intermediate value of α = 5 × 10⁻⁴, where α is the usual dimensionless diffusion coefficient, in units of c_s H, where c_s is the local midplane sound speed and H the disk's pressure scale height (Shakura & Sunyaev 1973). This value is consistent with observational data on turbulence in protoplanetary disks, suggesting α of the order of 10⁻³–10⁻⁴ (Pinte et al. 2016; Flaherty et al. 2017, 2018). The disk viscosity is a key parameter for the pebble problem, with smaller values of α leading to larger pebbles, thus generally faster inward drift, and longer persistence timescales of the gas locally enriched by pebble evaporation. All solid material is considered to be in the form of pebbles, with an initial particle size of 1 μm, which then evolve by growth and drift (e.g., Birnstiel et al. 2012).

The disk's resulting C/O values in the solid and gas phase are shown in the left panel of Figure 6. At t = 0 the disk C/O values reproduce our static setup. At larger times, however, the effect of drifting pebbles becomes noticeable very quickly. Pebbles drifting across the CO ice line will start enriching the gas phase in CO (also see Figure 7). Some of this gas diffuses outward again, condensing on the inward-drifting pebbles and increasing the pebble C/O value to unity just outside the CO ice line. The same effect is visible just outside the CO₂ and H₂O ice lines, where the solid C/O values reach 0.5 and 0, respectively. Away from the ice lines the C/O of the solids remains largely unchanged, however. At the same time we note that the solid surface density will drop significantly over the simulated time owing to pebble drift, by up to two orders of magnitude, while the gas surface density only drops by less than one order of magnitude. Inside the CO₂ ice line the C/O ratio of the gas immediately drops at t > 0 as CO₂ evaporates off the inward-drifting pebbles. At later times the gas's C/O starts rising again as the CO gas that has evaporated off the pebbles inside the CO ice lines reaches the inner disk regions, due to the disk's viscous evolution. An analogous evolution can be observed for the disk gas inside the H₂O ice line; the gas's C/O value first drops significantly, due to the water evaporating off the pebbles, but rises again at later times as gas enriched in CO₂ and CO viscously spreads inward.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For studying the effect of pebble accretion, drift, and evaporation, we investigated four scenarios with our formation inversion setup. Case (i) is simply applying the disk compositional model, as determined by the pebble drift and evaporation framework, in the formation forward model. Case (ii) is like Case (i), but putting an upper limit of 20 M_⊕ on the mass that can be accreted as solids, accounting for the concept of the pebble isolation mass (see Section 2.3 and Bitsch et al. 2018b). Case (iii) is like Case (ii), but replacing the upper mass limit on the planetary mass from the spectroscopic retrieval (M _P < 14 M_Jup; Mollière et al. 2020) with a tighter prior from the dynamical mass estimate reported in Brandt et al. (2021), that is, ${M}_{P}={9.6}_{-1.8}^{+1.9}\ {M}_{\mathrm{Jup}}$. Case (iv) is like Case (iii), but increasing the pebble isolation mass to 100 M_⊕. The reasoning for testing these different cases will be discussed below, where we summarize the inversion results obtained for the different cases.

The compositional fit for Case (i) is depicted by the leftmost ellipse in the bottom panel of Figure 4. In this scenario pebble drift, evaporation, and accretion are able to reproduce the observed abundance pattern of HR 8799e. Conceptually, Case (i) simply tests whether the results of the static disk model inversion change when introducing pebbles, but it does not yet apply any prior knowledge on how pebbles are accreted, such as the concept of the pebble isolation mass. The reason for the good compositional fit of Case (i) becomes evident when studying the right panel of Figure 6, which shows the resulting posterior for the most likely accretion location of the solids for HR 8799e. Because the solid C/O values do not change significantly, except for just outside the ice lines, the result that significant accretion of solids from outside the CO ice line is likely does not change when compared to the static setup of the disk composition. Just outside the CO ice line the probability goes down, however, because CO gas recondensing on the pebbles drives up the pebbles' C/O to values larger than the planetary one. What is noticeable is that the region inside the CO ice line at t > 0 is somewhat more likely when compared to t = 0. This is because the disk gas, enriched by CO from evaporating pebbles, and with C/O = 1, is of high enough metallicity to somewhat offset the C/O value of solids accreted inside the CO ice line, which is too low when compared to the planet. The enrichment of the disk gas in CO over time is shown in Figure 7. We note that the likelihood for accreting a significant amount of pebbles decreases over time because pebbles will drain to the inner parts of the disk. We neglect this effect here.

In Case (i) the upper limit on the planetary mass from the spectroscopic retrieval, together with a superstellar atmospheric metallicity, leads to a 1σ upper limit of 570 M_⊕ on the accreted pebble mass. Such a high value is inconsistent with the concept of the pebble isolation mass. Therefore, we deem Case (ii), where we set an upper limit of 20 M_⊕ on the accreted solid mass, a more likely scenario. We note that the pebble isolation mass is a function of the disk viscosity α and that it is very sensitive to the disk aspect ratio (M_iso ∝ [H/r]³, with H being the disk scale height). The value of 20 M_⊕ is what we derive for the HR 8799 disk model at the location of the CO ice line, using the scaling relations reported in Bitsch et al. (2018b). The compositional fit for Case (ii) is shown in the bottom panel of Figure 4 (second ellipse from the left). Also in this scenario pebble drift is able to reproduce the observed abundance pattern of HR 8799e but leads to a generally somewhat lower planetary metal enrichment. We note that these results assume that all accreted pebbles are visible in the atmosphere, which is equivalent to full core dissolution and mixing. Moreover, in order to allow pebble enrichment to have a noticeable effect on the planet composition, the inversion constrains the planetary mass to <3.5 M_Jup. What is more, as a result of the prior limit on the accreted solid mass and the planetary-mass prior, the inversion deems scenarios more likely where the composition of the accreted gas has more impact than in Case (i). The resulting probability distribution on the locations a_gas where the planet accreted its gas is shown in the left panel of Figure 8. The most likely locations and times for the gas accretion correspond to the situation where the gas enriched by the evaporating pebbles reaches approximately C/O values of 0.6, corresponding to the planet's atmosphere (see left panel of Figure 6). In this scenario the most likely formation accretion location is inside of HR 8799e's current orbital position, which would require some outward migration if taken at face value.

Zoom In Zoom Out Reset image size

In Case (iii) we study the effect of enforcing an upper limit on solid accretion, due to the pebble isolation mass, and a tighter constraint on the planetary mass. The mass prior stems from a dynamical analysis based on the orbital characterization of the HR 8799 system and accelerations from the Gaia-Hipparcos catalog (Brandt et al. 2021). The compositional fit for Case (iii) is shown in the bottom panel of Figure 4 (second ellipse from the right). In this case the inversion struggles to reproduce the observed enrichment pattern of HR 8799e; while it fits the atmospheric C/O ratio well, the planetary enrichment is generally too low, but it improves at later times. This is explained from the fact that the high mass prior assumed for HR 8799e, together with the low pebble isolation mass, does not allow for the pebbles to play a significant role in the planet enrichment, while especially at early times the disk gas is not enriched enough by gas that has evaporated off the inward-drifting pebbles. This situation is thus alleviated at later times, when the disk gas enrichment increases, but it is never enough to fully reproduce the planetary metal enrichment. The right panel of Figure 8 shows the probability distribution of a_gas. Because only gas accretion is able to affect the planetary composition noticeably in Case (iii), it is essentially a higher-contrast version of the a_gas distribution of Case (ii), shown in the left panel.

Case (iv) essentially studies the case when the planet started forming very far outside the CO ice line, in the outer parts of the disk. Due to the disk flaring, H/r increases toward the outer disk, and we would find M_iso = 50 M_⊕, corresponding to H/r = 0.07, at 200 au. Because we are interested in an upper limit on what pebble accretion could contribute, we also assume that the disk viscosity is very high for the M_iso calculation (α = 0.004, instead of the nominal 0.0005); this results in M_iso = 100 M_⊕. The corresponding enrichment pattern of the planet is shown in the bottom panel of Figure 4, rightmost ellipse. Unsurprisingly, not only the planet C/O but also its metal enrichment are better fit now, when compared to Case (iii), leading to a good fit overall. As expected, a_solid values outside the CO ice line are the most likely for this case, with some additional gas accretion from within the CO₂ ice line.

From our investigation it thus becomes evident that pebbles alone, for average M_iso values, may not be sufficient to fully explain the observed abundance pattern of HR 8799e, even when making the assumption that all pebbles accreted onto the planet (likely onto the forming planetary core) mix into the visible atmosphere. This conclusion hinges on at least three assumptions. First, if the planetary mass was actually lower than reported in Brandt et al. (2021), enriching the planet by the accreted solid pebbles becomes easier. This is shown by our Case (ii), where the inversion when imposing M_iso = 20 M_Jup resulted in a good fit by constraining the planetary mass to below 3.5 M_Jup. Next, if the pebble isolation mass is much higher than our baseline case (e.g., 100 instead of 20 M_⊕), which is possible for large disk viscosities and the planet initially forming far out in the disk, pebble enrichment becomes a likely scenario for explaining the abundance pattern of HR 8799 again. Lastly, if the composition of the disk is different from our baseline case, even M_iso = 20 M_⊕ with the dynamical mass prior of HR 8799e (Case iii) becomes a likely scenario again. This is seen in Figure 9, where we show what happens when running Case (iii) again, but assuming the λ Boo–type composition of HR 8799 for the disk composition. Due to the carbon and oxygen content being ∼30% higher in this case, the accretion of gas that is enriched by evaporated pebbles leads to a better agreement with the total atmospheric metal enrichment. We note that all of these conclusions are based on the spectroscopic retrieval result for HR 8799e, and a slightly lower retrieved metallicity would make pebble accretion more likely again. Due to the large uncertainties on the atmospheric metallicity, even the pebble scenario with the worst fit (Case iii) is only about one standard deviation away from the mean composition derived in the spectroscopic retrievals.

Zoom In Zoom Out Reset image size

Lastly, it should be kept in mind that other likely important effects connected to pebbles were not studied here. For example, we neglected the effects that the outer HR 8799 planets may have had on the pebble flux that reaches the inner disk and therefore HR 8799e's position. Outer planets may prevent pebbles from drifting inward and evaporating at the CO ice line (e.g., Bitsch et al. 2021; Schneider & Bitsch 2021a). It is unclear to what degree this effect is important here because the giant planets may have formed late enough that some pebble drift may already have taken place in the disk before shutting off the pebble flux. In addition, if HR 8799e, the innermost planet, formed first (high surface densities and orbital periods, and thus shorter accretion timescales), it may have been less affected by the formation of the outer planets.

Above it was studied how inferences drawn from a simple formation model change if chemical evolution of the protoplanetary disk, or the drift, evaporation and accretion of pebbles is included. As discussed in Section 2, planet formation is the combination of quite a number of key processes. A concurrent formation inversion with all the ingredients appears both numerically and conceptually unworkable, at the moment. It will still be instructive, however, to add certain aspects of the planet formation problem to such inversion calculations, to study their influence in isolation, or to assess the magnitude of their importance for atmospheric compositions. In Table 1 we list the way in which many of the aspects mentioned in Section 2 may be studied via inversion of the formation process.

To give an example, it would be straightforward to feed disk compositional models that include the disk's self-shadowing into the inversion framework. This process has been suggested by Ohno & Ueda (2021), where the shadowing is caused by a dust pileup at the water ice line. Depending on the grain properties and densities, such a scenario may allow for very volatile species such as CO, N₂, and even noble gases such as Ar to condense at distances from the star that are nominally too hot. To study such an effect, various midplane disk and abundance structures, for differing dust density contrasts, could be explored.

Another setup that would be instructive is to further investigate the effect of incomplete mixing between the deep interior and planetary atmosphere for gas giant planets. As discussed in Section 2, the metal enrichment inferred from atmospheric characterization studies is only a lower limit for the true planetary metal enrichment. Where available, an upper limit could be placed based on the analyses of planetary bulk metallicities, as obtained in Thorngren et al. (2016) and Thorngren & Fortney (2019). The impact of metallicity gradients could potentially be studied by adding a parameter f_mix that describes whether the metals accreted during formation fully mix (f_mix = 1) into the atmosphere or not (f_mix = 0). As long as a planetary atmosphere is of superstellar metallicity, f_mix will simply be inversely correlated with the accreted solid mass (if pebble evaporation is neglected). Once a planetary atmosphere is of stellar or substellar metallicity, f_mix may also correlate with the formation location of a planet, depending on the disk's abundance structure. An example for this can be constructed by considering our inversion results for HR 8799e in the static disk picture. Because the atmospheric metallicity is high and the planet has a stellar C/O value, the inferred atmospheric C/O ratio could only be reproduced by accreting solids from outside the CO ice line. If the planet's atmospheric metallicity was stellar, it could instead have formed at any location in the disk, as long as the location of gas accretion is equal to the location of solid accretion (measured with respect to the ice lines). If f_mix was added as a free parameter, small f_mix values would again have yielded regions outside the CO ice line as the most likely region of origin for the solids accreted by the planet.

It is also conceivable to construct a three-component model for the formation inversion that separates the planetary mass into three reservoirs: solids accreted onto the core, solids accreted and mixed into the gaseous envelope, and the gas itself. Each of these would also be associated with a parameter that described where the corresponding material was accreted. One could then define an f_mix for the core (or solids in the deep interior) that would describe the degree to which the deep core dissolves and mixes into the envelope.

Lastly, abundance constraints on the refractory content of a planet, as presented in Lothringer et al. (2021) for WASP-121b, for example, may be used to put an upper limit on the mass a planet accreted through pebbles. As long as the planet did not form very close around its star, where also refractories may enter the gas phase, refractories can only be incorporated into the planet by solid accretion. If the inferred amount of accreted refractories is higher than expected from the concept of the pebble isolation mass (see Section 2), an upper limit on the amount of accreted pebbles, as well as a lower limit on the amount of accreted planetesimals (or other impactors such as smaller planets; e.g., Ginzburg & Chiang 2020), may be constrainable. Similar constraints may be obtained from the cases where the volatile enrichment of a planet is higher even than what pebble evaporation in a disk may provide: any additional volatile mass must then be accreted in the form of volatile ices.

The importance of the planetary C/O ratio for informing planet formation has been discussed in Sections 1 and 2. In these sections we also discuss the complications that may make the picture likely more complex than suggested by the foundational study by Öberg et al. (2011). For example, while dominant solid enrichment is generally expected to lead to stellar or substellar C/O ratios and superstellar planetary metallicities, pebble drift and evaporation may lead to superstellar enrichments of the gas at C/O values both smaller or larger than stellar. Interestingly, however, superstellar C/O values and enrichments are difficult to obtain without considering pebble evaporation (also see Section 3.2), so a large population of planets with such abundance characteristics may indicate a dominant role of pebbles for setting planetary abundances. Similarly, a large enough overall metal enrichment of a planet, especially if formed in the outer disk, may be difficult to explain from pebble evaporation, even for low disk viscosities, which would point more toward planetesimal accretion playing an important role.

In general, C/O is also popular because it determines the relative abundances of the spectrally active C- and O-bearing molecules in the exoplanet atmospheres such as H₂O, CH₄, CO, CO₂, HCN, and C₂H₂. C/O therefore regulates the spectral appearance of a planet in the near- to mid-infrared (e.g., Fortney et al. 2005; Seager et al. 2005; Madhusudhan 2012; Moses et al. 2013; Mollière et al. 2015; Molaverdikhani et al. 2019a; Goyal et al. 2020; Hobbs et al. 2021b).

For reference, Figure 10 shows under which atmospheric conditions the absorbing species that trace the C/O ratio in gas-dominated planets may be visible. Also see, for example, Lodders & Fegley (2002) for a detailed description of the atmospheric chemistry. For temperatures below about 1000 K the atmosphere will be rich in H₂O and CH₄; for higher temperatures these species will be converted into CO until either C or O runs out, depending on the C/O ratio. For high temperatures and C/O ≳ 1, CH₄ will thus be visible; for high temperatures and C/O ≲ 1, H₂O will be visible. For further increasing temperatures and C/O ≳ 1, CH₄ is replaced by increasing amounts of C₂H₂ and $\mathrm{HCN}$ (see, e.g., Madhusudhan 2012; Mollière et al. 2015). We note that the chemical transitions mentioned here also depend on the local atmospheric pressure (Mollière et al. 2015; Molaverdikhani et al. 2019a). CO can still be visible in cool atmospheres, especially of self-luminous brown dwarfs and planets, because atmospheric mixing may transport CO-rich gas from the deep (hotter) atmosphere to the photosphere (e.g., Zahnle & Marley 2014; Miles et al. 2020). Whether such disequilibrium abundances are expected for irradiated (often transiting) planets is less clear because the insolation leads to more isothermal atmospheres. For planets that are still strongly cooling (with a high internal temperature) or heated by processes such as eccentricity dampening, CH₄ may be strongly suppressed (Fortney et al. 2020). As mentioned, ground-based high-contrast or high-resolution observations have started to obtain the first useful constraints on C/O. The state of the art will greatly improve once JWST and later ARIEL allow for a larger census of planetary compositions (also see our discussion in Section 1).

Zoom In Zoom Out Reset image size

The importance of atmospheric nitrogen-bearing species such as NH₃ or $\mathrm{HCN}$ for constraining exoplanet formation has been recognized recently, especially for planets forming in the outer solar and extrasolar disks. The reason for this is that nitrogen, predominantly in the form of N₂ in protoplanetary disks, is extremely volatile. Planets forming at increasingly larger distances, when dominated by solid metal enrichment, will therefore exhibit increasingly lower N/O or N/C ratios, and vice versa if dominated by gas metal enrichment. This is because several ice lines of C- and O-bearing species are crossed toward larger orbital radii, while N₂ stays in the gas phase (Turrini et al. 2021). If the planet forms at wide enough orbital distances, eventually N₂ will freeze out as well, leading to an enhanced atmospheric nitrogen content, which will scale similarly with atmospheric metallicity as the abundances of C- and O-bearing species. The high nitrogen content of Jupiter has therefore led to the interpretation that Jupiter formed in the outer regions of the solar system, beyond the location of the N₂ ice line at ∼30 au, which is also consistent with the planet's elevated abundance of noble gases (e.g., Owen & Encrenaz 2003; Bitsch et al. 2015; Bosman et al. 2019; Öberg & Wordsworth 2019; Cridland et al. 2020b). Similar to the discussion of C/O, the situation is likely more complicated also for the nitrogen enrichment. Both disk self-shadowing (see Ohno & Ueda 2021, and our discussion in Section 3.3) and pebble drift and evaporation (Schneider & Bitsch 2021b) are likely complicating factors. We also note that a planet that forms late within a protoplanetary disk's lifetime may be less sensitive to N₂, as cosmic-ray ionization may process N₂ to NH₃ ice over Myr timescales, such that the importance of NH₃ and its ice line increases over time, with the ice line of NH₃ being much closer to the star than that of N₂ (Semenov & Wiebe 2011).

In exoplanets the only spectrally active nitrogen-bearing species of relevance are NH₃ and HCN. N₂, which is the dominating nitrogen bearer at larger temperatures, has negligible opacity in the near- and mid-infrared. NH₃, on the other hand, should be detectable in the mid-IR using JWST in exoplanet atmospheres (e.g., Danielski et al. 2018). Moreover, evidence for NH₃ has been seen in high-resolution studies (Giacobbe et al. 2021; Sánchez-López et al. 2022). For C/O ≲ 1, NH₃ is only abundant up to ∼500 K (Lodders & Fegley 2002). HCN, on the other hand, will be visible for temperatures of 1500 K or larger, if C/O > 1 (e.g., Madhusudhan 2012; Mollière et al. 2015). We indicate these detectability ranges in Figure 10. Chemical disequilibrium may (or may not) allow for NH₃ or HCN to be visible at intermediate temperatures (500 K < T < 1500 K) in irradiated planets as well (see, e.g., MacDonald & Madhusudhan 2017; Fortney et al. 2020; Hobbs et al. 2021a, and references therein). For self-luminous planets disequilibrium chemistry may play less of a role for N, as isoabundance lines are parallel to atmospheric pressure–temperature profiles (Zahnle & Marley 2014). As before, all chemical transition temperatures also depend on the atmospheric pressure.

We note that the chemical behavior of N-, C-, and O-bearing species described here mostly hinges on chemical equilibrium or simple atmospheric disequilibrium treatments, also considering the planetary atmospheres to be one-dimensional and of mostly scaled solar abundances (except for the C/O ratio). The recent and intriguing results of Giacobbe et al. (2021), who detected H₂O, CO, $\mathrm{HCN}$, C₂H₂, NH₃, and CH₄ in the atmosphere of HD 209458b (with an equilibrium temperature of ∼1500 K), are a reminder that atmospheric chemistry may be much more complex than discussed above. An important effect is the horizontal advection of chemical abundances predicted from coupling chemical models to the output of 3D general circulation models (e.g., Agúndez et al. 2014; Baeyens et al. 2021), which could also be connected to condensate rainout (Sánchez-López et al. 2022). Photochemistry is also important, especially in the upper atmospheric layers (e.g., Kopparapu et al. 2012; Venot et al. 2012; Molaverdikhani et al. 2019b). Telescopes such as JWST, current high-resolution spectrographs, and ultimately instruments mounted on ELT-class telescopes will allow us to investigate these effects more thoroughly.

Measuring the refractory content of an atmosphere could provide unique insight into a planet's formation history. As has been argued recently by Lothringer et al. (2021), measuring the refractory-to-oxygen ratio R/O of a planet constrains the importance of metal enrichment by rocky accretion relative to icy or gaseous accretion. Here R stands for any element that traces the refractory content of the planet (Fe, Na, K, Si, Mg, Ti, ...), or an average of such elements. As argued further, this may allow the placement of constraints on whether the planet (or its solid building blocks) migrated significantly during formation. We argue that R/O may potentially even be useful to constrain the relative importance of pebble and planetesimal accretion, in the core accretion paradigm; see also Section 3.3 and the discussion in Schneider & Bitsch (2021b).

In Figure 10 we indicate the temperature ranges over which various refractory-tracing atmospheric absorbers are visible. We refer the reader to Appendix D for a discussion of the chemistry of the refractory-tracing absorbers. Lothringer et al. (2021) put emphasis on ultrahot Jupiters, for which various refractory elements exist as molecules (metal oxides or hydrides), atoms, or ions in the gas phase. We also note that species such as H₂S and PH₃ may be useful refractory tracers at intermediate atmospheric temperatures (Wang et al. 2017; Öberg & Wordsworth 2019). While H₂S and PH₃ are volatile species, the dominant carrier of P and S atoms in a protoplanetary disk appears to be refractory species (Öberg & Wordsworth 2019). Moreover, measuring the abundances of Na and K in planetary atmospheres may be worthwhile to trace of the refractory content (Welbanks et al. 2019).

Refractory cloud species may affect planetary spectra by muting molecular features and reddening the spectral energy distribution. Silicate particles like MgSiO₃ and Mg₂SiO₄ are especially interesting, as they may lead to visible absorption features around 10 μm (e.g., Cushing et al. 2006; Wakeford & Sing 2015). Due to the complex microphysical problem of cloud formation (e.g., Rossow 1978; Powell et al. 2018; Woitke et al. 2020), measuring a refractory abundance from observed cloud absorption may prove difficult, however. Moreover, clouds may complicate measuring and interpreting the abundances of gas refractory species due to cold trapping by condensate rainout (e.g., Spiegel et al. 2009; Parmentier et al. 2016). An interesting alternative to silicate clouds could be searching for the absorption of gaseous SiO at ∼7 μm with JWST's MIRI instrument. SiO is promising because it is the most abundant Si-bearing gas species after the silicates evaporate (Visscher et al. 2010) and is more stable than H₂O against dissociation (by about 500 K). This should allow detection of this species in ultrahot Jupiters.

Other metal oxides such as TiO, VO, AlO, and CaO have features in the optical and near-infrared (e.g., Sharp & Burrows 2007; Gandhi & Madhusudhan 2019; Lothringer et al. 2020). Similarly, metal hydrides may be useful refractory tracers, at similar temperatures to the metal oxides. Species such as FeH, CaH, MgH, NaH, CrH, and TiH all have absorption features in the optical and near-infrared (Sharp & Burrows 2007; Gandhi & Madhusudhan 2019). Metal atoms are visible in the atmosphere once the refractory clouds are no longer present (e.g., Mg, Fe), or once the dominant molecular species (such as SiO for Si) have been dissociated. Mg, Fe, Ca, Cr, Ni, V, Na, and maybe Co have been detected in the ultrahot atmospheres of KELT-9b and WASP-121b in the optical (Hoeijmakers et al. 2019, 2020). Finally, metal ions become visible in the hottest atmospheres as soon as the atoms have been ionized. This has led to the detection of Fe+, Ti+, Cr+, Sc+, Y+, and maybe Sr+ in the hottest known exoplanet KELT-9b in the optical (Hoeijmakers et al. 2019).

Here we describe the disk setup used for inverting the Öberg et al. (2011) model (see Section 3.1.1), which assumes a static disk composition. More specifically, our disk model is based on the description in Öberg & Wordsworth (2019), which assumes a static power law for the disk temperature and density of the young solar nebula. Inside its ice line position a given volatile species is in the gas phase, whereas outside it is in the solid phase. For every species the mass fraction compared to the total disk mass is tabulated, in addition to the mass fractions of the constituent atoms within a volatile species. We also account for refractory material, which we include in our framework by setting its ice line position to zero. The background gases H₂ and He are included by setting their ice line to 1000 au. This ensures that the refractory and background species stay condensed/gaseous within the simulation domain. Table 2 lists the mass fractions, ice line positions, and atomic composition of all considered disk species and their constituent atoms. We note that the ice line positions given here are those expected for the disk around HR 8799, which have been obtained from our ANDES disk model; see Section 3.2.

The mass fractions were obtained using the provisional protosolar nebula composition from Öberg & Wordsworth (2019) with some modifications; see below. The abundance of a species i is given as number fractions x_i in Öberg & Wordsworth (2019), compared to the number of hydrogen atoms n_H. This was converted into mass fractions m_i relative to the total disk mass using

$$\begin{eqnarray}&&{m}_{i}=\displaystyle \frac{{\mu }_{i}{x}_{i}}{1.4},\end{eqnarray} \tag{ B1 }$$

with μ_i being the molecular mass of species i in atomic mass units. This expression is obtained from the disk's total mass M_tot, from the mass of species M_i of species i, and by setting

$$\begin{eqnarray}&&{m}_{i}=\displaystyle \frac{{M}_{i}}{{M}_{\mathrm{tot}}}\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&\approx \,\displaystyle \frac{{\mu }_{i}{x}_{i}{n}_{{\rm{H}}}}{{\mu }_{{\rm{H}}}{n}_{{\rm{H}}}+{\mu }_{\mathrm{He}}{n}_{\mathrm{He}}}\end{eqnarray} \tag{ B3 }$$

$$\begin{eqnarray}&&=\,\displaystyle \frac{{\mu }_{i}{x}_{i}}{{\mu }_{{\rm{H}}}+{\mu }_{\mathrm{He}}{n}_{\mathrm{He}}/{n}_{{\rm{H}}}},\end{eqnarray} \tag{ B4 }$$

where it was assumed that most of the disk mass is contributed by H and He atoms. Setting μ_H = 1 and μ_He = 4 and assuming that n_He/n_H = 0.1 (see Table 8 in Lodders 2019, for the recommended protosolar abundances) leads to the relation given in Equation (B1).

The refractory composition model was likewise constructed using the information given in Öberg & Wordsworth (2019). In their model, this results in 30% of all oxygen in the form of refractory silicates, identical to the amount of oxygen in H₂O. We assumed that the silicates consist purely of MgSiO₃, which also conserves the solar Mg/Si abundance ratio, which is close to unity (Asplund et al. 2009; Lodders 2019). The refractory carbon component plus volatile organics account for 50% of all carbon atoms, with a 3:1 ratio between the two. To simplify, we added the carbon of the volatile organics by increasing the CO abundance, taking the required oxygen from the water mass reservoir. Because the fate of organic carbon, especially in the inner part of the disk, is uncertain anyway (see, e.g., Mordasini et al. 2016; Cridland et al. 2019, and references therein), we decided to forgo a more careful treatment of the organic carbon reservoir for our conceptual study here. Iron, sulfur, and phosphorous atoms were assumed to only be present in the refractory phase. The C, Fe, S, and P abundances, relative to H, were taken from Asplund et al. (2009). Our resulting C/O ratio distribution for the disk solid and gas components is shown in Figure 1.

Here we describe the ANDES chemistry model used for inverting the formation model including disk chemical evolution (see Section 3.1.2), which includes the evolution of the disk's chemical composition. In ANDES, surface reactions are described by the Langmuir–Hinshelwood mechanism and are not limited to hydrogenation. H and H₂ tunneling through reaction barriers is also included. Any dynamical effects on the distribution of C and O, such as drifting grains, are omitted. The surface density profile is described by a power law with the exponent equal to 1.5. The chemical network is based on ALCHEMIC (Semenov et al. 2010; Semenov & Wiebe 2011), with updated binding energies from Cuppen et al. (2017). It incorporates the effects of extreme ultraviolet irradiation, cosmic rays, and radionuclides as ionization sources. The dust size distribution is described by a power law with p = −3.5 between 0.005 and 25 μm, which reflects dust growth in disks compared to the interstellar medium. It is used to calculate the radiation field and dust temperature in the disk's upper layers. An average grain radius of 0.35 μm is used for calculating surface reaction rates. The disk abundances are initialized assuming that all volatile species are in the gas phase, using the same abundances as used for the static disk model; see Table 2.

H₂S condenses into NH₄SH at ∼200 K; the higher-temperature condensates MnS, ZnS, and Na₂S are of minor importance (Lodders 2010). For C/O < 1, H₂S dissociates at ∼2000 K, while it moves into species like CS for C/O > 1 at temperatures around 1500–2000 K.

PH₃ condenses into H₃PO₄ at 500 K. Its presence in Jupiter's atmosphere at lower temperature indicates a deep quenching point, however, such that PH₃ may still be visible at lower temperatures (see, e.g., Baudino et al. 2017). For temperatures approaching 1000 K, PH₃ is increasingly converted into PH₂.

Na condenses into Na₂S at ∼900 K. In principle, alkalis such as Na could also be sequestered into high-temperature condensates such as feldspars. However, this likely does not occur owing to the rainout of silicates (e.g., Line et al. 2017) that deplete Si from the atmosphere, which is needed for feldspar formation. Above 900 K Na is thus in the gas phase, until it gets ionized at around 2500 K.

K condenses into KCl at ∼900 K. In analogy to Na, sequestration of K into feldspars likely does not occur owing to silicate rainout. Thus, K is in the gas phase from ∼900 to ∼2000 K, after which it is ionized. Ionization occurs at temperatures roughly 500 K cooler than for Na.

As mentioned above, the refractory cloud species Na₂S and KCl likely form at temperatures below 900 K. Here we focus on the remaining cloud species forming at intermediate to hot temperatures and concentrate on those carrying the largest mass and/or opacity, often using the data given in Wakeford et al. (2017), or our own equilibrium chemistry calculations. Refractory clouds can exist if the atmospheric temperature is below their respective evaporation temperature. Like all chemical transitions discussed here, this temperature depends on the elemental abundance and local atmospheric pressure. Under our adopted standard conditions silicates such as MgSiO₃ and Mg₂SiO₄ evaporate at temperatures around ∼1600 K, while iron clouds are stable until ∼1700 K. VO and calcium titanates are stable until 1600–1800 K, respectively. Aluminum-bearing condensates such as Al₂O₃, which are among the most stable ones, evaporate around 1900 K. Among the species listed here, potentially only silicates, Al₂O₃, and KCl may actually form in the visible part of the atmosphere, as these species have low surface energies, leading to high nucleation rates (Gao et al. 2020). The cloud bases will reside deeper inside the planetary atmosphere for lower temperatures, with the cloud particles entering from above, or settling below the photosphere. For brown dwarfs this temperature-dependent removal of silicate clouds is thought to cause the L–T transition, which typically occurs at T_eff = 1200–1400 K (e.g., Best et al. 2021, and references therein), while for planets and low-gravity brown dwarfs this limiting temperature may be as low as approximately 1000 K (e.g., Marley et al. 2012; Morley et al. 2012; Charnay et al. 2018).

SiO is an especially interesting molecule for tracing the abundance of the refractory silicates in the atmosphere such as MgSiO₃ or Mg₂SiO₄. As soon as the silicates evaporate (around 1600 K), their constituent atoms move into the gas phase. While atomic Mg is then the preferred gaseous form of Mg, Si will move into SiO (Visscher et al. 2010). For C/O ≳ 1, SiO enters the gas phase at around 1300 K, which is when SiC evaporates. SiO then starts to dissociate around 3500 K for C/O ≲ 1, while moving into SiS around 2000 K for C/O ≳ 1. The local evaporation-dependent temperatures given here could be lower than the observed transition as a function of planetary effective temperature, where high-pressure cloud formation could cold-trap Si into silicates.

Similar to SiO, the other metal oxides form as soon as the refractory clouds evaporate. Possible species of interest are TiO, VO, SiO, AlO, and CaO (e.g., Sharp & Burrows 2007; Gandhi & Madhusudhan 2019), with the relevant evaporation temperatures of the clouds ranging from ∼1600 to 1900 K at our adopted standard conditions. Again, these are then only expected to be visible in the atmosphere if not cold-trapped into condensates at lower altitudes, that is, higher pressures. Except for SiO (see discussion in the SiO section above), most of these metal oxides are not expected to form in atmospheres with C/O ≳ 1 (Madhusudhan 2012; Gandhi & Madhusudhan 2019). In general, metal oxides will be destroyed by dissociation at high enough temperatures, with TiO and VO dissociating at temperatures similar to water (around 3000 K). As stated above, SiO is a bit more stable, dissociating from temperatures higher by about 500 K.

Similar to metal oxides, metal hydrides such as FeH, CaH, MgH, NaH, CrH, and TiH may form as soon as the refractory clouds have been evaporated, thus at local atmospheric temperatures of around 1600–1900 K and cooler temperatures for NaH, as Na₂S evaporates at ∼900 K already. Of course, the cold trapping statement from above holds here as well. The hydrides will be destroyed by thermal dissociation at high enough temperatures, for example, MgH and FeH dissociate around 3000 K or so (Lothringer et al. 2018). For these two species it should also be noted that the main gas-phase carrier is atomic Mg and Fe in (ultra)hot Jupiter atmospheres and that MgH and FeH are less abundant by about a factor of 10⁴ (Visscher et al. 2010).

Metal atoms can exist in the gas phase as soon as the sequestering refractory condensates evaporate (modulo cold trapping). As mentioned above, gaseous Fe and Mg are the main gas carriers of these elements once the silicates and iron condensates are gone. Si takes over as the main gas species only after SiO is dissociated. Fe, Mg, Ti, Ca, and Ni all are ionized between 3500 and 4000 K or so, with Al being ionized at somewhat lower temperatures (Kitzmann et al. 2018; Lothringer et al. 2018; Hoeijmakers et al. 2019).

Metal ions become visible in the atmosphere as soon as the atoms have been ionized; see immediately above.